ofsp2026 20_simple | Aerosand.cn

Created time

Mar 28, 2026 08:36 AM

type

Post

status

Published

date

Mar 28, 2026

slug

ofsp2026 20_simple

summary

0. 前言

我们回忆计算流体力学的基本控制方程 https://aerosand.cc/docs/cfd/cfdb/05_general-conservation/

质量方程

动量方程

参考 03_momentumConservation，可以简化成通用形式

能量方程

总结有通用形式基本方程

在 OpenFOAM 中，SIMPLE（Semi-Implicit Method for Pressure Linked Equations）算法用于求解稳态问题。

本文主要讨论

密度处理和压力参考

压力速度耦合方程

非线性项的显隐

SIMPLE 算法

SIMPLE 代码框架

1. 密度和压力

为了方便讨论，考虑无重力稳态不可压缩流动的 NS 方程

连续方程（质量方程）

通用形式的动量方程（包括一个出于简单考虑暂时为零的源项 S ）

在 OpenFOAM 的不可压缩求解器中，压力 p 实际上是动压头，隐含的已经除以了密度，即实际上

其中，P 是物理压力，此时的 p 的单位为

即使已经使用了动压头，动压头也是基于参考点的相对压力。因为对于流体动力学的计算来说，压力梯度才会真正影响流动。而设置真实压力，可能会因为数量级悬殊会引起计算精度的损失。一般会选择压力参考。参考点一般选择在计算域内部（例如 cell 0），尽量保证计算域内相对压力的计算稳定性。

类似的，在 OpenFOAM 不可压缩求解器中，粘性力项也默认除以了密度，用运动粘度 nu 代替了动力粘度 mu ，即

2. 控制方程

控制方程有

连续方程（质量方程）

动量方程为

对于此控制方程，在求解中面临两个问题：

的非线性

压力和速度的耦合

我们需要一定的算法来求解此方程组。

3. 非线性

回忆之前的对流项，物理量 \phi 和质量通量 F 作为整体被称为对流通量。而这里相当于是自己输运自己（速度输运）。

📌

Note

如果完整考虑密度的话，其实应该称为“动量输运”。

现在的对流项，作体积积分后离散

对于其中的非线性部分 (UU)，实践中更倾向于线性处理，也就是令质量通量中的速度，使用上一时间步的值作为已知值，即所谓的显式处理。质量通量之外的速度，用待求时间步的值作为未知量，即所谓的隐式处理。当然，这样会导致计算过程中，非线性项存在一个延迟。不过这对于稳态计算来说并不重要。

对于瞬态计算来说，仍然可以忽略非线性延迟而作一个线性化处理，也可以单独对非线性项作迭代进行处理。但是，单独作迭代会增加计算消耗。当时间步很小的时候，连续解之间的差别就会很小，此时非线性延迟也就不再显著。

讨论可以知道，对流项中的质量通量是变化的核心，可以作为一个独立的量进行处理。

📌

Note

需要注意，既然进行了线性化处理，可以想到上式中 a_{P},a_{N} 对某一时间步的计算来说可以作为已知量参与计算。但是对于全部时间步变化来说，这两个系数是跟着速度更新而变化的。这些系数也构成了下文中的 M 矩阵（所以M矩阵在每个时间步都是不同的）。

实践中，一般使用 SIMPLE 压力-速度耦合算法来计算稳态流动计算，使用 PISO 压力-速度耦合算法来计算瞬态流动计算。后续将分别介绍这些算法。

4. 压力速度耦合

动量方程的形式可以简化为（忽略其他源项）

可以想见场的系数矩阵 M 是一个对角占优的稀疏方阵（或者说我们希望它能够对角占优），M 矩阵的对角线都是离散方程中的本元素，同一行上，对角线前后（在该行也许不紧挨）的元素都是和该单元相邻的单元。

4.1. 动量预测

最基本的思路是，在某个时间步或者初始时间步，我们由上一步求得的已知速度压力场或者初始已知速度压力场，根据动量方程直接求出一个预测的速度场。

约定每个迭代步【动量预测】求解动量方程得到的预测速度表示为

动量方程为

为了方便理解，我们摘抄历史版本的主要代码如下 https://github.com/OpenFOAM/OpenFOAM-2.2.x/blob/master/applications/solvers/incompressible/simpleFoam/UEqn.H

这个直接求解动量方程的步骤被称为【动量预测】，求得预测速度 U^{pre}。

📌

Question

为什么压力构建不放入 UEqn 的构建内呢？

读者可以先自行思考，我们后续将详细讨论。

📌

Note

也许会有读者疑惑，整体上看速度和压力同样都是要求解的。为什么这里构建的 UEqn 求解的一定是速度而不是压力呢？

实际上 fvm:: 离散项返回了一个矩阵，即构建成了 fvVectorMatrix 类型的一个矩阵，离散操作对象就是参数 U，也就是待求未知量。而 fvc:: 离散项返回了一个场（用于计算的已知量），即计算出了 volVectorField 类型的一个场。

4.2. 压力修正

基于【动量预测】得到的预测速度，同时需要满足质量方程的约束，可以以此修正压力。

我们先处理系数矩阵 M 。

📌

Note

记得 M 矩阵是由类似于上文讨论的离散后中的系数构成的。

从 M 矩阵中分解出对角矩阵 A （不是求解线性代数系统 Ax=b 中的 A），对角矩阵 A 可以很容易求出逆矩阵。

OpenFOAM 中的对角矩阵以及求解方法已经高度抽象化，也更易读

对角矩阵 A 的逆矩阵在代码中为

动量方程的左侧抽离对角矩阵，可以写成

相应的非对角矩阵为

OpenFOAM 中也提供可以直接调用的成员方法 UEqn.H() 。

我们回到动量方程。分解后的动量方程为

两边同时乘以 A 的逆矩阵

其中的也被称为 HbyA(U)（其实就是 H 除以 A 的意思），在代码中为

继续数学推导有

速度还需要满足连续方程

即有

整理有【压力修正方程】

理论上，为了求解得精确压力，我们应该提供精确的。

实际上，我们只能提供基于预测速度的参与求解计算。

这样的操作，在实际上是假设忽略对计算的影响并不大。

📌

这样一来，通过【动量预测】得到的预测速度来计算新的压力（修正压力）